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Observations from old forests underestimate 
climate change effects on tree mortality 
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Understanding climate change-associated tree mortality is central to linking climate change 
impacts and forest structure and function. However, whether temporal increases in tree 
mortality are attributed to climate change or stand developmental processes remains 
uncertain. Furthermore, interpreting the climate change-associated tree mortality estimated 
from old forests for regional forests rests on an un-tested assumption that the effects of 
climate change are the same for young and old forests. Here we disentangle the effects 
of climate change and stand developmental processes on tree mortality. We show that both 
climate change and forest development processes influence temporal mortality increases, 
climate change-associated increases are significantly higher in young than old forests, and 
higher increases in younger forests are a result of their higher sensitivity to regional warming 
and drought. We anticipate our analysis to be a starting point for more comprehensive 
examinations of how forest ecosystems might respond to climate change. 
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Observational studies have relied on old forests to quantify 
temporal trends in background tree mortality, and have 
shown that tree mortality has increased with recent global 
warming, increasing atmospheric C0 2 and decreasing water 
availability in tropical, temperate and boreal forests 1-5 . Old 
forests are assumed to be in an equilibrium state where tree 
mortality is matched by recruitment, and temporal changes in 
tree mortality at the forest level due to endogenous processes are 
assumed to be weak at this stage of stand development 6 ' 7 . 
However, whether there is ever equilibrium is a contentious 
issue 8 . Although the effects of competition estimated by temporal 
changes of stand basal area (SB A) and stand density on tree 
mortality have been investigated in some previous studies 4,5 , 
conspecific negative density dependence and tree aging could also 
influence tree dynamics 9 ' 10 . Furthermore, the temporal increase 
of tree mortality could be driven by stand development processes, 
but not necessarily climate change 11,12 . Unfortunately, previous 
studies have used unsuitable statistical methods that marginalize 
either climate or non-climate drivers for longitudinal data in 
which these drivers are highly correlated 13 . Therefore, the relative 
roles that climate change and stand development processes have 
on the temporal changes in tree mortality remain uncertain. 

Furthermore, climate change- associated tree mortality 
increases in old forests have been used to represent the tree 
mortality response of regional forests to climate change 14 , based 
on the assumption that effects of climate change on tree mortality 
are same for young and old forests. However, this assumption has 
not been specifically tested. Verification of the assumption is 
essential for disturbance-driven ecosystems such as boreal forests. 
Because of high fire frequency, the boreal forests are a mosaic of 
stands at various developmental stages with old forests 
accounting for only a small portion of the landscape 15 ' 16 . With 
a predicted increase of fire frequency associated with global 
warming, the mean age of the boreal forests is expected to 
decrease, as will the portion of old forests 17 . Thus, understanding 
how young forests respond to climate change is essential to 
predict how boreal forests will respond to future climate change. 
If tree mortality responses to climate change differ between young 
and old forests, increased tree mortality based on old forests will 
be biased for regional forest predictions. 

Here we analyzed tree mortality patterns of boreal tree species 
in western Canada using 887 permanent sample plots measured 
between 1958 and 2007. To account for uncertainties in sampling, 
models and parameters, we used Bayesian models 18, to 
disentangle endogenous (stand development) and exogenous 
(year, temperature anomaly or drought) effects on individual tree 
mortality for five major boreal tree species, Populus tremuloides 
Michx., Populus balsamifera L., Pinus banksiana Lamb., Picea 
mariana Mill, and Picea glauca (Moench) Voss. We tested the 
assumption that tree mortality increases associated with climate 



change do not differ with forest age. We also examined whether 
tree mortality increases are linked to recent regional warming and 
its negative consequence on water availability. We show that both 
climate change and forest development processes influence 
temporal mortality increases, climate change-associated 
increases are significantly higher in young than old forests, and 
higher increases in younger forests are a result of their higher 
sensitivity to regional warming and drought. Furthermore, 
climate change-associated mortality increases are stronger for 
moist-habitat- adapted Populus balsamifera and late-successional 
Picea mariana and Picea glauca than Populus tremuloides and 
Pinus banksiana. 



Results 

Effects of year and stand development processes on tree mor- 
tality. With all data pooled, the effects of year, used to represent 
climate change drivers as a whole in both Model 1 (with only year 
as a predictor) and Model 2 (with year, endogenous processes and 
their interactions as predictors) were significantly positive for all 
study species (Table 1; Fig. 1), indicating significant increases of 
annual mortality probability during the study period (Fig. 2). The 
two models, however, produced different estimates for climate 
change-associated tree mortality: year effects estimated by Model 
1 are bigger for Populus tremuloides, Pinus banksiana, Picea 
mariana and Picea glauca, but smaller for Populus balsamifera 
than those by Model 2 (Fig. la), resulting in different estimates of 
annual tree mortality probability from the two models (Fig. 2). 

Tree mortality was strongly affected by stand development 
processes and their interactions with year (Fig. lb, Supplementary 
Table S4). Year effect decreased significantly with stand age (SA) 
for all study species (Table 1). Year effect also increased with 
stand crowding as measured by SBA and in stands with more 
conspecific individuals for Picea mariana and Picea glauca 
(Table 1). Tree mortality increases associated with year differed 
among species with stronger year effects on late-successional 
Picea mariana and Picea glauca and moist-habitat- adapted 
Populus balsamifera than Populus tremuloides and Pinus bank- 
siana (Fig. la). 

When data were analyzed separately for young forests (initial 
SA < 80 years) and old forests (initial SA > 80 years, same as in 
Peng et a/. 5 ), year effect was consistently higher in young forests 
than old forests for all species (Fig. la), supporting the strong 
declining year effect with increasing SA (Table 1). Predicted 
annual mortality increased by 2.88%, 4.70%, 5.64% and 6.18% per 
year in young forests, but increased at the lower rates of 1.65%, 
2.24%, 3.89% and 4.12% in old forests for Populus tremuloides, 
Populus balsamifera, Picea mariana and Picea glauca, respectively 
(Supplementary Fig. S3). For drought- tolerant Pinus banksiana, 
annual mortality increased by 2.09% per year in young forests, 



Table 1 | Effect of year and its interactive effects with stand development processes on annual mortality probability. 



Model Term 



Populus tremuloides 



Populus balsamifera 



Pinus banksiana 



Picea mariana 



Picea glauca 



Model 1 Intercept 

Year ( x10" 2 ) 
AUC 

Model 2 Intercept 

Year ( x10" 2 ) 
Year x SA ( x10 
Year x 0.5 RBA 
Year x SBA ( x10" 
Year x rFSBA ( x 10 
AUC 



-4.03 (-4.10 to -3.98) 
2.63 (2.46-2.80) 
0.711 

-4.33 (-4.40 to -4.27) 
2.42 (2.13-2.71) 
4 ) -4.67 (-5.54 to -3.82) 
NS 

" 4 ) NS 
NS 
0.812 



-3.54 (-3.63 to -3.44) 
2.64 (2.34-2.95) 
0.748 

-3.67 (-3.81 to -3.59) 

3.63 (3.24-4.01) 
-4.60 (-6.06 to -3.13) 
NS 
NS 
NS 
0.819 



-5.18 (-5.37 to -5.00) 
3.41 (2.59-4.23) 
0.768 

-5.44 (-5.63 to -5.25) 

1.85 (0.80-2.89) 
-3.40 (-4.63 to -2.17) 
NS 
NS 
NS 
0.825 



-4.94 (-5.08 to -4.80) 
4.83 (4.39-5.27) 
0.750 

-4.76 (-4.90 to -4.62) 

3.83 (3.14-4.49) 
-3.16 (-4.41 to -1.94) 
NS 

5.46 (1.02-9.97) 
3.95 (2.20-5.69) 
0.777 



-5.06 (-5.17 to -4.95) 
4.47 (4.26-4.67) 
0.754 

-4.86 (-4.97 to -4.75) 

3.93 (3.55-4.33) 
-5.10 (-5.85 to -4.35) 
NS 

9.32 (6.85-11.70) 
6.05 (5.14-7.01) 
0.776 



Abbreviations: AUC, area under the receiver operating characteristic curve; NS, corresponding predictor's posterior 95% credible interval covers 0 in the full model and the predictor was removed in the 
reduced model (Methods); RBA, relative basal area; rFSBA, ratio of focal species' basal area to SBA; SA, stand age; SBA, stand basal area. 

Values are estimated parameters (mean and 95% credible interval in brackets). Full table for Model 2 is presented in Supplementary Table S4. The models predictive performances were evaluated by 
AUC. 
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Figure 1 | Year effect on annual tree mortality probability and sensitivity scores of predictors, (a) Year effect on annual tree mortality probability, 
logit (p), estimated by Model 1 (without endogenous factors as predictors) and Model 2 (with endogenous factors as predictors). Models were 
separately developed all plots (All), young plots (Young, initial SA <80 years) and old plots (Old, initial SA >80 years), respectively. Error bars are 
95% credible intervals, (b) Sensitivity scores. For each species and age group (All, Young or Old), sensitivity scores of predictors from Model 1 are 
on the left and Model 2 on the right. 



Populus balsamifera 
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1960 1970 1980 1990 2000 1960 1970 1980 1990 2000 1960 1970 1980 1990 2000 

Year 

Figure 2 | Predicted temporal trends of annual mortality probability associated with calendar year. The predicted means (solid lines) and their 
95% credible intervals (dotted lines) of annual mortality probability are derived by using equation exp(/0 — 1, in which is the fitted year coefficient from 
Model 1 (red) and Model 2 (blue) for each respective species in Table 1. 



but year had no detectable effect on its mortality in old forests 
(Supplementary Fig. S3). In both young and old forests, 
endogenous factors were critical drivers of tree mortality for all 
study species (Fig. lb; Supplementary Tables S6 and S7). 



mortality of all species and greater effects for Populus balsamifera, 
Picea mariana and Picea glauca than the other two species. 
Furthermore, the ACMIA effects were stronger in young than old 
forests for all species except Pinus banksiana (Table 2). 



Warming and drought effects on tree mortality. For the study 
area, annual temperature anomaly (ATA) increased at 
0.034 °C year _ and annual climate moisture index anomaly 
(ACMIA) decreased at 0.088 year - l , whereas the growing season 
precipitation anomaly slightly increased at 0.761 mm year ~ 1 
between 1958 and 2007 (Supplementary Fig. S4; Supplementary 
Table S7). To examine whether increase of mortality could be 
attributed to regional warming and drought, we developed 
models that replaced year by these climatic variables (Models 3 
and 4). The ATA models showed positive main ATA effects on 
mortality of all species and higher ATA effects in young than 
old forests of all species except Pinus banksiana (Table 2). 
The ACMIA models indicated negative main ACMIA effects on 



Discussion 

Unlike previous attributions of temporal increases in tree 
mortality to either climate change 1- or stand development 
processes 10 ' 12 , by using Bayesian models, we show that both stand 
development processes and climate change have affected 
temporal increases in tree mortality. These findings are not 
likely a result of methodological problems because heterogeneity 
in sampling strategies had no effect on mortality (Supplementary 
Figs S5 and S6), though we note that we have not explicitly 
considered alternate exogenous factors here. Tree mortality 
associated with year decreased for all species except Populus 
balsamifera when stand development factors were included as 
predictors. The differences in year effect between the models 
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Table 2 | Effect of ATA or A CM I A and its interactive effects with stand development processes on annual mortality probability. 



Model Term 



Populus tremuloides 



Populus balsamifera 



Pinus banksiana 



Picea mariana 



Picea glauca 



Model 3 Intercept 
ATA ( 
ATA 
ATA 
ATA 
ATA 
AUC 
Intercept 
ACMIA ( 
ACMIA : 
ACMIA : 
ACMIA : 
ACMIA : 
AUC 



10" 1 ) 

SA (x10" 3 ; 
0.5 RBA 
SBA ( x 10 ~ 
rFSBA(x10" 



Model 4 



<10" 2 ) 
SA ( x 10~ 3 ) 
0.5 RBA (x10" 
SBA 

rFSBA (xlCr 



-4.29 (-4.36 to -4.23) 

4.67 (4.19-5.15) 
-9.26 (-11.06 to -7.53) 

NS 

NS 

1.88 (0.27-3.47) 
0.801 

-4.30 (-4.37 to -4.23) 
-3.59 (-4.54 to -2.64) 

1.93 (1.56-2.29) 
-8.63 (-13.09 to -4.25) 
NS 
NS 
0.815 



-3.91 (-4.03 to -3.79) 

5.17 (4.35-6.00) 
-12.80 (-15.84 to -9.77) 
NS 
NS 
NS 
0.743 

-3.83 (-3.95 to -3.70) 
-7.27 (-8.84 to -5.70) 
1.99 (1.39-2.60) 

NS 

NS 

7.33 (1.79-12.85) 
0.798 



-5.44 (-5.64 to -5.24) 
5.44 (3.35-7.55) 
NS 

-1.55 (-2.49 to -0.59) 
NS 
NS 
0.774 

-5.45 (-5.66 to -5.24) 
-1.58 (-0.14 to -3.02) 

NS 

NS 

NS 

NS 
0.803 



-5.00 (-5.11 to -4.89) 

3.97 (3.00-4.90) 
-4.37 (-6.57 to -2.33) 
NS 
NS 
NS 
0.745 

- 4.61 ( - 4.79 to - 4.44) 
-7.77 (-9.72 to -5.83) 
0.57 (0.12-1.02) 

NS 

NS 

-12.01 (-17.81 to -6.23) 
0.759 



-4.90 (-4.98 to -4.81) 

5.26 (4.71-5.81) 
-7.52 (-8.91 to -6.17) 
NS 

7.75 (2.00-13.52) 
9.46 (7.39-11.49) 
0.728 

_4.68 (-4.82 to -4.56) 
-4.44 (-5.37 to -3.48) 
1.17 (0.89-1.46) 

NS 

NS 

-14.57 (-18.41 to -10.74) 
0.783 



Abbreviations: ACMIA, annual climate moisture index anomaly; ATA, annual temperature anomaly; AUC, area under the receiver operating characteristic curve; NS, corresponding predictor's posterior 
95% credible interval covers 0 in the full model and the predictor was removed in the reduced model (Methods); RBA, relative basal area; rFSBA, ratio of focal species' basal area to SBA; SA, stand age; 
SBA, stand basal area. 

Values are estimated parameters (mean and 95% credible interval in brackets). Full tables for Model 3 and Model 4 are presented in Supplementary Tables S8 and S9, respectively. The models predictive 
performances were evaluated by AUC. 



provide support for the notion that studies excluding endogenous 
factors can potentially produce biased estimates of climate change 
effects in boreal forests 13 . The decreased year effect on tree 
mortality is attributable to stand development processes during 
the study period. For example, both SA and SBA increased 
(Supplementary Table S10), both of which positively affected tree 
mortality (Supplementary Table S4). Also, an increase in the 
conspecific density of shade-tolerant Picea mariana and Picea 
glauca (Supplementary Table S10), reflecting the nature of 
secondary succession in boreal forests 6 ' 20 , intensified a negative 
density- dependence effect on these two species (Table 1). The 
increased year effect in moist-habitat-adapted Populus 
balsamifera suggests that reduced water availability might have 
overweighed endogenous effects. 

Our results show that mortality increases associated with 
climate change were significantly higher in young than old 
forests, and the higher increase in younger forests appears to be 
due to higher mortality sensitivity to recent regional warming and 
the resulting negative consequence on water availability. This 
finding indicates that climate change- associated increases in tree 
mortality would be underestimated if only old forests are used to 
represent regional forests. Compared with old forests, even-aged 
young forests established after a stand- replacing disturbance may 
experience greater competition for space and nutrients among 
young trees that tend to occupy the same ecological niche 6, , 
resulting in them being more vulnerable to external stressors such 
as climate change-associated drought 21 ' 22 . The positive effect of 
regional warming on mortality could also be attributed to 
the interdependencies among warming, drought and forest 
pests 23-25 . In our study area, outbreaks of forest pests could 
have a role in the temporal mortality increases 26-28 . 

Mortality responses to climate change differed among species, 
suggesting that the regional forest may be undergoing forest 
compositional shift that is independent of endogenous forest 
succession. The highest mortality increase in Populus balsamifera 
among pioneer species indicates that reduced water availability by 
regional warming has the greatest influence on species adapted to 
moist habitats. The higher mortality increases of late-successional 
Picea mariana and Picea glauca than those of early- successional 
Populus tremuloides and Pinus banksiana in both young and old 
forests (Figs 1 and 2) suggest that the regional forest will likely 
become further dominated by early- successional species if current 
warming trends continue. The different responses of tree 
mortality among species, coupled with those by climate change- 
associated fire activities 15-17 ' that promote early- successional 
species 30 , may pose a significant conservation challenge for the 
region under further global warming. 

4 



We show evidence that long-term tree mortality trends are 
influenced by both stand development processes and long-term 
climate trends, that is, regional warming and drought. Note that 
the strength of these drivers on tree mortality differs among 
species, and their influences interact both among endogenous 
factors and between endogenous and exogenous factors. Because 
of the higher sensitivity to recent climate change in young than 
old forests, climate change-associated tree mortality could be 
underestimated if mortality estimates from old forests are used to 
represent regional forests. Tree mortality responses to climate 
change differed among species, and these different responses can 
have strong implications for future forest composition in the 
western boreal forest. 

Methods 

Study area and the forest inventory data. We analyzed longitudinal data 
from 887 permanent sampling plots established in Alberta and Saskatchewan 
(Supplementary Data 1 and Supplementary Fig. SI). Details for plot establishment, 
measurements and selection criteria are described in Supplementary Methods. 

Explanatory variables. To examine whether long-term tree mortality trends are 
affected by long-term climate change trends, similar to previous studies 1-5 ' 12 , we 
used the middle calendar year of a census period, during which tree mortality 
measurements were made between two successive censuses, to represent climate 
change drivers as a whole. To account for the effects of endogenous factors 
including asymmetric competition, stand crowding, tree ageing and inter- specific 
interaction on tree mortality, we used the tree's relative basal area (RBA; ratio of 
subject tree basal area to the mean tree basal area of the stand), SBA, SA and 
the ratio of focal species' basal area to SBA, calculated using the preceding 
measurement of each census period 10 . 

Climate anomalies were denned as the departure from the long-term climate 
means 31 . Because our dependent variable, tree status, was observed at the end of a 
census period, the departure from the long-term mean for each climate variable 
was denned as the difference between its mean value during the census period 
and its long-term mean. Long-term climate means were calculated based on the 
1950-2007 time period, during which plot measurements were taken. Three sets of 
climate anomalies were calculated: ATA, growing season precipitation anomaly 
and ACMIA. We derived temperature and precipitation data using scale-free 
ClimateWNA 32 ' 33 and climate moisture index (CMI) using BioSIM 34 from 
latitude, longitude and elevation of each sample plot. Annual CMI was the sum of 
monthly CMIs over 12-month periods from last 1 August to 31 July of the current 
year, where monthly CMI was the quantity of monthly precipitation minus 
monthly potential evapotranspiration, which was computed using a simplified 
form of the Penman-Monteith equation 35 . A smaller CMI value indicated a drier 
condition. The summary statistics for explanatory variables were presented in 
Supplementary Table S2. 

Annual mortality probability calculations. We used individual tree mortality 
analyses rather than plot-level mortality analyses to accommodate the individual 
tree-level variable, RBA, because (1) RBA is the strongest predictor of tree mortality 
in non-equilibrium boreal forests and (2) its control on tree mortality can change 
with SBA and ratio of focal species' basal area to SBA 10 . Similar to previous 
studies 36 ' 37 , we were interested in estimating annual tree mortality probability 
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where measurement intervals varied (Supplementary Fig. S2e). We derived annual 
mortality probability from tree status as follows. 

The survival probability, S, of tree i, for census period j of t years, in plot k is 
expressed as a function of annual survival probability, s: 

S* = >t (!) 

Annual mortality probability, pyk = 1 — Sy& (2) 

The mortality probability for a census period, = 1 — S,^ (3) 

Lastly, we related a tree's status at the end of the census period (My*, alive or 
dead) to annual mortality probability using Bernoulli likelihood: 

My fc ~ Bernoulli (1 - (1 - p ijk ) fi *) (4) 

Statistical analyses. For each species, we developed individual tree mortality 
models using Hierarchical Bayesian logistic regression with plot identity as a 
random effect to account for uncertainties potentially affected by different forest 
characteristics or in different climatic and/or edaphic conditions 18 ' 19 . 

logit(Py fc ) = a + PxX ijk + n k n ~ N(0, o\) (5) 

where a and P are the intercept and estimated coefficients, respectively; Xy k represents 
explanatory variables corresponding to tree i for census period j in plot k; and the 
model included a term n k to describe the random effect of sampling plots. As a rule of 
thumb, n k is in normal distribution with a mean of 0. 

To examine whether inclusion of endogenous factors alters climate change- 
associated tree mortality, we compared the year effect between models with and 
without inclusion of endogenous factors as predictors, that is, Model 1 and 
Model 2, respectively. 

logit(pyfc) = a + P x Yearyjt + n k (Model 1 ) 

logit(py*) = a + f(En) + p u x Yeary fc + 0 12 x Year^ x 0.5 RBA * 

+ p 13 x Yeary fc x SBA yfc (Model 2) 

+ p iA x Yeary* x SA^ + ft 5 x Year yfc x rFSBA^ + n k 

where fiEn) in Model 2 was developed based on our previous work 10 to account for 
the endogenous effects on tree mortality: 

f(En) = logit(py fc ) = p i x0.5 mA >* + P 2 xSBA ijk 

+ P 3 x SAy, + £ 4 x rFSBAy, + £ 5 x 0.5**** x SBA yfc 

+ p 6 x 0.5^ x SAy fc + p 7 x 0.5**** x rFSBAy* 

+ £ 8 x SBAy fc x SAy fc + P 9 x SBAyjt x rFSBAy* + j& 10 x SA^ x rFSBAy fc 

where fts are the coefficients. We modelled the mortality probability as an exponential 
function of RBA because our previous results showed that mortality probability 
decreases exponentially with RBA 10 . The Akaike information criterion also indicated 
that the exponential transformation of RBA gave better fits for all species 
(Supplementary Table S3). Furthermore, our previous work 10 showed near-linear 
mortality trends with tree aging and with SBA after taking account of competition and 
species interaction. Consequently, we modelled mortality probability as a linear 
function of SA and SBA. 

Both Models 1 and 2 were developed for all plots, young forest plots and old 
forest plots, respectively, to examine how the inclusion of endogenous factors affect 
year effect in these three scenarios. Following Peng et al. 5 , we denned that young 
forests were < 80 years of age and old forests > 80 years of age at the first 
census. 

To examine whether increases of mortality could be attributed to regional warming 
and its negative consequence on water availability, we developed models that replaced 
year with the climatic variables in Model 2, resulting in Models 3 and 4. 

Iogit(py0 = a + f(En) + p n xATAy* + P l2 x ATA ijk x0.5 RBA * 

+ P l3 x AT Ay* x SBAy* + P u X ATA ijk x SA ijk (Model 3) 

+ P l5 x ArA ijk x rFSBA y * + n k 

logit(py fc ) = a +f(En) + p u xACMIAy* + p l2 x ACMIAy* x0.5 RBA <* 
+ P 13 x ACMIAy- fc x SBA^ 

+ p u x ACMIAy fc x SAy fc + p i5 x ACMIAy fc x rFSBAy fc + n k 

(Model 4) 

To examine the temporal trends of endogenous and exogenous factors at the plot 
level, Y, we used the Hierarchical Bayesian linear (HBL) model. 

Y jk = a + PxYearjk + n k n ~ N(0, o 2 k ) (7) 

where Yj k is dependent variable, that is, endogenous factors and exogenous factors; a 
and P are intercept and estimated coefficients, respectively; Year ;fc represents middle 
calendar year at for census period j in plot k; and the model included a term n k to 
describe the random effect of sampling plots. 



For all analyses, the Bayesian Markov Chain Monte Carlo methods were 
implemented using /AGS 38 called from R 39 with rjags package 40 . All coefficients were 
assigned non- informative priors. All independent variables were centred to reduce 
their correlations and speed up convergence. For each model, we evaluated 
convergence by running two independent chains with different initial values and 
monitoring the Gelman-Rubin statistic 41 . When convergence was confirmed, an 
additional 10,000 iterations with thinning of half were used to calculate the mean, s.d. 
and 95% credible interval for each coefficient from the posterior distribution. 

We developed full models and sequentially removed dimensionless explanatory 
variables whose posterior 95% credible interval covered 0. The reduced models, 
assessed by Deviance Information Criterion, were better than or similar to full models 
(Supplementary Table Sll). Consequently, we selected the reduced models as our 
final models. To assess the adequacy of the final models, we calculated the area under 
the receiver operating characteristic curve (AUC) based on the mean of each 
coefficient using ROCR package 42 . A value of AUC > 0.8 indicates that a model has 
excellent discriminatory power, and a value > 0.7 indicates good discriminatory 
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power . 

Finally, we calculated the mortality sensitivity scores to endogenous factors, 
exogenous factors and their interactions as a measure of the influences of these factors 
on mortality. Similar to Dietze and Moorcroft 37 , we denned sensitivity score on 
annual mortality probability as the s.d. of the predicted annual mortality probability 
for each of the above three groups of factors, holding the other two groups of variables 
at their mean. 

To present how annual mortality probability changes over the study 
period, we summarized posterior distribution of annual mortality probability 
(mean and 95% credible interval) using mean of each parameter (/?) from the final 
models. We also calculated the annual fractional change of mortality probability 
associated with year over the study period using the equation exp(P) — 1, 
in which P is fitted coefficient for Year. 

Assessing possible methodological problems. We first examined the spatial 
auto -correlation and heterogeneity in sampling strategies following the method 
described by van Mantgem et al. . For all study species, Mantel tests showed no 
evidence of spatial auto -correlation among the plots for annual mortality 
probability, annual fractional change of mortality probability or the effect of plot 
identity on annual mortality probability (Supplementary Table S12). The selected 
permanent sampling plots were located in ecologically heterogeneous regions, and 
close geographic configuration does not imply similarity in forest characteristics. 
Our result that mortality probability increased more so in young forests than old 
forests could also be mistaken if the young forest plots were spatially clustered. 
Thus, we also examined the spatial auto -correlation for initial SAs. For all study 
species, the mantel tests showed no evidence of spatial auto -correlation among the 
plots for initial SA (Supplementary Table S12). 

The linear regressions that related annual mortality fractional change to plot size 
and average census interval showed that neither plot size nor average census 
interval was significantly related to annual mortality fractional change 
(Supplementary Figs S5 and S6). These results indicate that the variation of 
mortality changes could not be attributed to plot size heterogeneity or the variation 
in census intervals. 
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